clc;
clear;

animate = false; % 是否绘制动画 true or false

a = 0; b = 2; % 区域 [a,b]
T_0 = 0; % 初始时间
x_int = 0.5; % 波前的初始位置
alpha = 1; % 参数
c0 = x_int - alpha*T_0; % (x_int,T_0) 特征线
h = 0.01; % 空间网格步长
lambda = 0.5; % 网格比 tau/h
dt = lambda*h; % 时间网格步长
N_t = ceil(1/dt); % 时间网格剖分数
x = a:h:b; % 空间网格点坐标
N_x = length(x); % 空间网格剖分数

% 精确解
u_e = zeros(N_x, N_t);
for n = 1:N_t
    for j = 1:N_x
        if x(j) <= c0 + alpha*n*dt
            u_e(j,n) = 1;
        else
            u_e(j,n) = 0;
        end
    end
end

% 初始化矩阵，用于存储各种数值解
u_FTFX = zeros(N_x, N_t);
u_FTBX = zeros(N_x, N_t);
u_FTCX = zeros(N_x, N_t);
u_LaxFriedrichs = zeros(N_x, N_t);
u_LaxWendroff = zeros(N_x, N_t);
u_BeamWarming = zeros(N_x, N_t);
u_leapFrog = zeros(N_x, N_t);

% 初值 u_0
for j = 1:N_x
    if x(j) <= x_int
        u_FTFX(j,1) = 1;
        u_FTBX(j,1) = 1;
        u_FTCX(j,1) = 1;
        u_LaxFriedrichs(j,1) = 1;
        u_LaxWendroff(j,1) = 1;
        u_BeamWarming(j,1) = 1;
        u_leapFrog(j,1) = 1;
        u_leapFrog(j,2) = 1;
    else
        u_FTFX(j,1) = 0;
        u_FTBX(j,1) = 0;
        u_FTCX(j,1) = 0;
        u_LaxFriedrichs(j,1) = 0;
        u_LaxWendroff(j,1) = 0;
        u_BeamWarming(j,1) = 0;
        u_leapFrog(j,1) = 0;
        u_leapFrog(j,2) = 0;
    end
end

% 边值 u(a,t)=1, u(b,t)=0
for n = 1:N_t
    % 左边界
    u_FTFX(1,n) = 1;
    u_FTBX(1,n) = 1;
    u_FTCX(1,n) = 1;
    u_LaxWendroff(1,n) = 1;
    u_LaxFriedrichs(1,n) = 1;
    u_BeamWarming(1,n) = 1;
    u_BeamWarming(2,n) = 1; % u(j-2,n)项
    u_leapFrog(1,n) = 1;
    % 右边界
    u_FTFX(N_x,n) = 0;
    u_FTBX(N_x,n) = 0;
    u_FTCX(N_x,n) = 0;
    u_LaxWendroff(N_x,n) = 0;
    u_LaxFriedrichs(N_x,n) = 0;
    u_BeamWarming(N_x,n) = 0;
    u_leapFrog(N_x,n) = 0;
end

% 计算各个格式的数值解
for n = 1:N_t
    for j = 1:N_x-1
        u_FTFX(j,n+1) = u_FTFX(j,n) - alpha*lambda*(u_FTFX(j+1,n)-u_FTFX(j,n));
    end
    for j = 2:N_x
        u_FTBX(j,n+1) = u_FTBX(j,n) - alpha*lambda*(u_FTBX(j,n)-u_FTBX(j-1,n));
    end
    for j = 2:N_x-1
        u_FTCX(j,n+1) = u_FTCX(j,n) - 0.5*alpha*lambda*(u_FTCX(j+1,n)-u_FTCX(j-1,n));
        u_LaxFriedrichs(j,n+1) = 0.5*(u_LaxFriedrichs(j+1,n)+u_LaxFriedrichs(j-1,n)) - 0.5*alpha*lambda*(u_LaxFriedrichs(j+1,n)-u_LaxFriedrichs(j-1,n));
        u_LaxWendroff(j,n+1) = u_LaxWendroff(j,n) - 0.5*alpha*lambda*(u_LaxWendroff(j+1,n)-u_LaxWendroff(j-1,n)) + 0.5*alpha^2*lambda^2*(u_LaxWendroff(j+1,n)-2*u_LaxWendroff(j,n)+u_LaxWendroff(j-1,n));
    end
    for j = 3:N_x
        u_BeamWarming(j,n+1) = u_BeamWarming(j,n) - alpha*lambda*(u_BeamWarming(j,n)-u_BeamWarming(j-1,n)) - 0.5*alpha*lambda*(1-alpha*lambda)*(u_BeamWarming(j,n)-2*u_BeamWarming(j-1,n)+u_BeamWarming(j-2,n));
    end
end
for n = 2:N_t
    for j = 2:N_x-1
        u_leapFrog(j,n+1) = u_leapFrog(j,n-1) - alpha*lambda*(u_leapFrog(j+1,n)-u_leapFrog(j-1,n));
    end
end


if ~animate % 不绘制动画
    figure;
    hold on;
    xmin = min(x);
    xmax = max(x);
    ymin = min(u_e(:,1));
    ymax = max(u_e(:,1));
    pad = 0.05;
    xmin = xmin - pad*(xmax-xmin);
    xmax = xmax + pad*(xmax-xmin);
    ymin = ymin - pad*(ymax-ymin);
    ymax = ymax + pad*(ymax-ymin);
    axis([xmin xmax ymin ymax]);
    grid on;
    n = ceil(N_t/2);
    title(['Time: ', num2str(n*dt)]);
    xlabel('x');
    ylabel('u(x,t)');
    plot(x, u_e(:,n), '-r');
    plot(x, u_FTFX(:,n), '-b*');
    plot(x, u_FTBX(:,n), '-g');
    plot(x, u_FTCX(:,n), '-c*');
    plot(x, u_LaxFriedrichs(:,n), '-m');
    plot(x, u_LaxWendroff(:,n), '-b');
    plot(x, u_BeamWarming(:,n), '-k');
    plot(x, u_leapFrog(:,n), '-y');
    legend({'Exact', 'FTFX', 'FTBX', 'FTCX', 'Lax-Friedrichs', 'Lax-Wendroff', 'Beam-Warming', 'Leap Frog'});
else % 绘制动画，保存视频文件
    videoObj = VideoWriter('FDM_stability_1D.mp4', 'MPEG-4');
    open(videoObj);
    figure;
    hold on;
    xmin = min(x);
    xmax = max(x);
    ymin = min(u_e(:,1));
    ymax = max(u_e(:,1));
    pad = 0.05;
    xmin = xmin - pad*(xmax-xmin);
    xmax = xmax + pad*(xmax-xmin);
    ymin = ymin - pad*(ymax-ymin);
    ymax = ymax + pad*(ymax-ymin);
    axis([xmin xmax ymin ymax]);
    grid on;
    for n = 1:N_t
        title(['Time: ', num2str(n*dt)]);
        xlabel('x');
        ylabel('u(x,t)');
        plot(x, u_e(:,n), '-r');
        % plot(x, u_FTFX(:,n), '-b*');
        plot(x, u_FTBX(:,n), '-g');
        % plot(x, u_FTCX(:,n), '-c*');
        plot(x, u_LaxFriedrichs(:,n), '-m');
        plot(x, u_LaxWendroff(:,n), '-b');
        plot(x, u_BeamWarming(:,n), '-k');
        plot(x, u_leapFrog(:,n), '-y');
        % legend({'Exact', 'FTFX', 'FTBX', 'FTCX', 'Lax-Friedrichs', 'Lax-Wendroff', 'Beam-Warming', 'Leap Frog'});
        legend({'Exact', 'FTBX', 'Lax-Friedrichs', 'Lax-Wendroff', 'Beam-Warming', 'Leap Frog'});
        frame = getframe(gcf);
        writeVideo(videoObj, frame);
        clf;
        hold on;
        axis([xmin xmax ymin ymax]);
        grid on;
    end
    close(videoObj);
    close(gcf);
end
